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Abstract 

A recently developed efficient recursive approach for analytically calculating the short-time evolution of the 
one-particle propagator to extremely high orders is applied here for numerically studying the thcrmody- 
namical and dynamical properties of a rotating ideal Bose gas of ^''Rb atoms in an anharmonic trap. At 
first, the one-particle energy spectrum of the system is obtained by diagonalizing the discretized short-time 
propagator. Using this, many-boson properties such as the condensation temperature, the ground-state 
occupancy, density profiles, and time-of-flight absorption pictures are calculated for varying rotation fre- 
quencies. The obtained results improve previous semiclassical calculations, in particular for smaller particle 
numbers. Furthermore, we find that typical time scales for a free expansion are increased by an order of 
magnitude for the delicate regime of both critical and overcritical rotation. 

Key words: Bose-Einstein condensation. Effective action. Exact diagonalization, Time-of-fiight 
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1. Introduction 

Bose-Einstein condensation (BEC) represents a 
macroscopic quantum phenomenon of recent broad 
research interest [Ji HH, Since its first experi- 
mental realization in 1995, it has been extensively 
studied experimentally, analytically, and numeri- 
cally. The two main research directions are weakly- 
interacting dilute gases in magneto-optical traps 
and strongly-interacting quantum gases in optical 
lattices. The behavior of a Bose-Einstein conden- 
sate under rotation is essential for understanding 
many fundamental BEC phenomena. For instance, 
its response to rotation represents one of the sem- 
inal hallmarks of superfiuidity [5|. However, once 
a harmonically trapped Bose-Einstein condensate 
is rotated critically, i.e. the rotation frequency be- 
comes so large that it compensates the radially con- 
fining harmonic frequency, the system turns out to 
be radially no longer trapped. In the absence of ad- 
ditional potential terms the condensate would start 
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to expand perpendicular to the rotation axis. For 
an overcritical rotation, this expansion would even 
be accelerated by the presence of a residual centrifu- 
gal force. In order to reach experimentally this del- 
icate regime of critical or overcritical rotation. Fet- 
ter suggested in Ref. Q to add an additional quartic 
term to the harmonic trap potential. Using a Gaus- 
sian laser beam propagating in the z-direction, this 
has been realized experimentally in Paris by Dal- 
ibard and co-workers for a BEC of 3 • 10^ atoms of 
®''Rb 0, Q. The resulting axially-symmetric trap 
with a small quartic anharmonicity in the xy-plane, 
seen by individual atoms, has the form 



-^-^Rb / 2 o2\ 2 , ^^Rb 2 2,^4 

Vbec = j'-i + ^-w^z +-r^ 



, (1) 



with the perpendicular radius rj_ — \/ +2/^, as 
well as the trap frequencies = 27r x 64.8 Hz, 
iOz = 27r X 11.0 Hz, and the trap anharmonic- 
ity fc = fcBEC = 2.6 X 10"" Jm-''. Further- 
more, the rotation frequency $7, which is measured 
in units of i.e. it is expressed by the ratio 
r = fi/wj^, represents the tunable control param- 
eter which could be experimentally increased up 
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to r = 1.04. This highest possible rotation fre- 
quency seems to coincide with an instabihty which 
fohows from a Thomas-Fermi sohition of the Gross- 
Pitaevskii equation [§]. 

As long as we can ignore the presence of two- 
particle interactions and approximately describe 
the system with the ideal Bose gas, its many- 
particle properties in the grand-canonical ensemble 
are exclusively derivable from one-particle states. 
When considering the thermodynamic limit, usu- 
ally the semiclassical approximation is applied, 
where the one-particle ground state Eq is retained 
and treated quantum mechanically, while all one- 
particles states above Eq are approximately treated 
as a continuum fid\. This semiclassical approxi- 
mation remains reasonable good irrespective of the 
rotation frequency Q once the total particle num- 
ber N is large enough and the trap anharmonic- 
ity k small enough. The latter condition implies 
that the underlying one-particle potential ^ has 
a small curvature around its minimum. However, 
in this context the question arises for which sys- 
tem parameters such a semiclassical approximation 
is not sufficient for a precise description of BEG 
phenomena, as well as when it finally breaks down, 
requiring a full quantum mechanical treatment of 
the system. 

In order to analyze this fundamental problem 
more quantitatively, it is mandatory to deter- 
mine the one-particle energy eigenvalues and eigen- 
functions fully quantum mechanically. To this 
end we apply a recently developed ultra-fast con- 
verging path-integral approach for calculating the 
imaginary-time evolution amplitude of gene ral non- 
relativistic many-body systems 11, 12, 13, 3, 
In this approach a hierarchy of discretized effec- 
tive actions is introduced where higher terms, com- 
pared to the naive action, substantially reduce er- 
rors in calculations of transition amplitudes. In par- 
ticular, this allows a systematic derivation of dis- 
cretized effective actions of level p, which lead to a 
1/iVP-convergence of discretized time-sliced transi- 
tion amplitudes within the continuum limit iV — > cx) 
of infinitely-many time slices. In addition, the im- 
proved convergence of transition amplitudes, calcu- 
lated with the level p effective action, can be used to 
construct higher-order analytic approximations for 
short-time propagators. Furthermore, the = 1 
time-slice approximation turns out to be valid for 
short imaginary times and useful for small values of 
the inverse temperature in applications in quantum 
statistical physics. The recent Refs. 16|,ll7[ demon- 



strate that this path-integral approach allows to 
determine a huge number of both eigenvalues and 
eigenfunctions of quantum systems with very high 
accuracy using exact numerical diagonalization. 

In this Letter we show how this path-integral 
approach is applied for studying both global and 
local properties of fast-rotating Bose-Einstein con- 
densates. To this end we proceed as follows: Sec. [5] 
briefly reviews the main ingredients of the path- 
integral approach in the general context of ideal 
Bose-Einstein condensates. Then we calculate in 
Sec. [3] a large number of energy eigenvalues and 
eigenfunctions for the anharmonic one-particle po- 
tential (dl). Afterwards, Sec. U discusses how a finite 
number of numerically available energy eigenvalues 
affects the results and how they can be improved 
by introducing systematic semiclassical corrections. 
On the basis of this precise numerical one-particle 
information. Sec. [5] studies global properties of a 
rotating condensate, for instance the condensation 
temperature Tc as a function of rotation frequency 

and the ground-state occupancy Nq /N as a func- 
tion of temperature T. Finally, Sec. [B] is devoted to 
the calculation of local properties of the condensate, 
such as density profiles and time-of-fiight absorp- 
tion pictures. Sec. [7] briefly summarizes the main 
results presented in this Letter. 

2. Many- Versus One-Particle Physics 

At first we demonstrate that a precise numerical 
access to one-particle eigenstates allows the calcu- 
lation of the condensation temperature and other 
thermodynamic properties of an ideal Bose gas. 
Afterwards, we briefiy review details of the path- 
integral effective action approach for a numerical 
study of a ^^Rb BEG in the anharmonic trap ([T]). 

2.1. Ideal Bose Gas 

For high temperatures, the grand-canonical par- 
tition function of an ideal Bose gas is given by 



(2) 



where v enumerates all possible configurations of 
the system, /3 — l/ksT represents the inverse tem- 
perature, and /Li denotes the chemical potential. As 
the ideal bosons do not interact, the system energy 
El, can be expressed in terms of single-particle en- 
ergy eigenvalues 



E,, 



7 , Af^(n) En , 



(3) 



where n counts single-particle energy states, while 
Nv{n) = 0, 1, 2, . . . and En stand for the occupancy 
and the energy eigenvalue of level n, respectively. 
Correspondingly, the number of particles in the sys- 
tem reads 

(4) 



Thus, the grand-canonical free energy J- 
— (lnZ)//3 results to be 



-/3(£;„-m) 



(5) 



where a subsequent Taylor expansion of the loga- 
rithm yields 



1 ^ e™^^ 

^ m=l 



m 



(6) 



usually denoted as the cumulant expansion. Thus, 
the many-body thermodynamic potential ^ of an 
ideal Bose gas is exclusively determined by single- 
particle states via the one-particle partition func- 
tion: 

Zi(/3) = Ve-'^^". (7) 



In principle, the above outlined exact calculation 
of the many-body free energy allows a full numeri- 
cal description of ideal Bose gases, and can be also 
applied for studies of dilute Bose gases in the case 
when interactions are negligible. However, it be- 
comes numerically very involved even for simple 
trapping potentials at low temperatures. In addi- 
tion to this, the BEC phase transition is achieved 
only in the thermodynamic limit of an infinite num- 
ber of atoms, thus making numerical studies of the 
condensation increasingly difficult. Usually, this 
problem is solved by fixing the chemical potential 
H at the low temperatures of the condensate phase 
to the ground-state energy, i.e. by setting = Eq, 
and to treat the ground state separately, by explic- 
itly taking into account its macroscopic occupation 
A^o- Thus, for low enough temperatures the grand- 
canonical free energy ^ is modified to 



T = 



1 -A e™^^ 
8 ^ m 

^ m=l 

+No{Eo''fi). 



-m/3-Bol 



(8) 



In order to avoid any double-counting, we have 
subtracted in the first line the contribution of the 



ground state within the one-particle partition func- 
tion, whereas the second line takes into account 
a possible macroscopic occupation of the ground 
state. The resulting total number of particles N = 
—dF/dfi follows to be 



N 



OO 



'"^^ [Zi{ml3) 



(9) 



This particle number equation serves different pur- 
poses in the respective phases. Within the gas 
phase, where the macroscopic occupation of the 
ground state vanishes, i.e. we have A^o = 0, 
Eq. ^ determines the temperature dependence of 
the chemical potential fj,. On the other hand, within 
the BEC phase the chemical potential /x coincides 
with its minimal value, i.e. the ground-state energy 
Eq, so Eq. © yields the temperature dependence 
of A^o- Therefore, the value of f3c — l/fcs7c, which 
characterizes the boundary between both phases, 
follows from Eq. © by setting Nq = and jjb = Eq: 



m— 1 



1] 



(10) 



We conclude that, for a given number N of ideal 
bosons, the condensation temperature can be ex- 
actly calculated only if both the single-particle 
ground-state energy Eq and the full temperature 
dependence of the one-particle partition function 
([7]) are known. 

2.2. Path-Integral Approach 

The recently developed path-integral approach 
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converging numerical calculation of all one-particle 

[ll,[l3|. Note that 



properties of quantum systems 
this general numerical approach is suitable to treat 
arbitrarily-shaped trap potentials and is also appli- 
cable for exact studies of many-body problems. In 
this Letter we consider BEC experiments performed 
by Jean Dalibard's group 0, Q , where a trapping 
potential with a quartic anharmonicity was real- 
ized. For this reason, we will specialize our numer- 
ical approach in the following to such potentials in 
order to quantitatively analyze this seminal exper- 
iment. 

The quantum statistical imaginary-time transi- 
tion amplitude A{a.,h;t) = (b|e~*^/''|a) follows in 
the continuum limit oo of infinitely-many time 
slices 

v4(a,b;t)= lim ^Ar(a,b;i). (11) 



Here the time-sliced amplitude AN{a,h;t) is a 
product of short-time amplitudes corresponding to 
the introduced time steps, and is expressed as a 
multiple integral of the function g-SN/ti^ where Sjy 
is usually called discretized (Euclidean) action [18| . 
For the simplest one-particle theory with the Eu- 
clidean Lagrangian L = q^/2 + V{q), where the 
coordinates are rescaled so that the mass is equal 
to unity, the naive discretized action is given by 



S 



N 



N-1 

E 



(12) 



with the abbreviations e = t/N, Sn — Qn+i — Qn, 
Qn = (Qn+i + qn)/2 and boundary conditions qo = 
a, q^v = b. Using the naive discretized action leads 
to the convergence of the transition amplitudes to 
the continuum as slow as 1/iV. 

Starting from the one-particle Schrodinger equa- 
tion in d spatial dimensions, we have derived a hi- 
erarchy of discretized effective actions S"^-* with a 
systematically improved convergence to the contin- 
uum result. The general short-time amplitude is 
first written in the form 



A(q„,q„+i;£) 



(27rfe)'i/2 



F-tH'(q„,5„;e) 



(13) 



were M^(q„,(5„;e) stands for the ideal discretized 
effective action, giving exact transition amplitudes 
in any discretization. In fact, the above equation 
is valid for any propagation time t, not just for 
short times, and if we were able to calculate ex- 
actly the ideal effective potential W, we could use it 
directly to calculate the long-time transition ampli- 
tude A{a, b; t). This is practically possible only for 
a limited number of exactly solvable models. How- 
ever, we have derived a recursive approach which al- 
lows for an efficient analytical calculation of a high- 
order short-time expansion of the effective potential 
[15|. If we insert the expansion T^(P-i) calculated 
to order e^^^ into Eq. p^ . we get the short-time 
amplitude 



A(P)(q„,q„+i;e) 



(27rfe)'^/2 



which is designated by p due to the fact that T/t^(P~i) 
is multiplied by an additional factor of e in the ex- 
ponent. This yields a result for the exponent which 



is correct up to order e^, i.e. its error is propor- 
tional to eP~^^. If we take into account the pre- factor 
l/(27rfe)'^/2, the above expression for the short time 
amplitude A'-p') gives overall errors proportional to 
gP-i-i/2 in d = 1. Consequently, in d = 2 spatial 
dimensions, which is the case relevant here, the er- 
rors of short-time transition amplitudes would be 
of the order e^. However, if we consider the N 
dependence of the above expression, the errors are 
always proportional to l/N^^^, since the pre-factor 
l/(27rfe)''/2 will be absorbed in the normalization 
of the amplitude A(a, b; t). When used to calculate 
a long-time transition amplitude, the product of N 
short-time amplitudes will subsequently lead to er- 
rors proportional to l/N^j as demonstrated conclu- 
sively earlier in Refs. [lil, [11 [H, El . 

Typical values of the inverse temperature /? in 
BEC experiments are quite small compared to the 
typical energy scale which is defined by the har- 
monic trap frequencies. For example, in the Paris 
experiment JJ the dimensionless value of h/3uj± 
ranges between 10~'^ and 10~^. Therefore, one 
can immediately use the above formula for the am- 
plitude A'-p'> and calculate the corresponding one- 
particle partition function by numerically integrat- 
ing the diagonal amplitude ^(a, a, /3) over the co- 
ordinate a. For small enough /? the above formula 
converges rapidly, and the amplitudes can be cal- 
culated exactly for all practical purposes. We refer 
to this approximation as the = 1 approximation, 
since it corresponds to having only one time-slice in 
the standard definition of the transition amplitude 
in the path-integral formalism. 

However, the direct use of this approach has sev- 
eral disadvantages. First of all, one still has to per- 
form an integral over the diagonal coordinate a in 
order to calculate the partition function. Second, 
this has to be done repeatedly for each value of the 
inverse temperature /?. And most importantly, one 
has also to extract the value of the ground-state 
energy in view of the particle number equation (|9]). 
In principle, this is done by studying the high-/? 
regime, where the short-time expansion (jl4p is not 
valid. Although this procedure works also for lower 
values of /3 |19j . it requires the numerical calcu- 
lation of the one-particle partition function and a 
detailed study of its dependence on the inverse tem- 
perature in order to obtain the ground-state energy 
with sufhcient precision. For this reason, the al- 
gorithm becomes numerically complex and difficult 
to use, especially in cases where the ground state is 
degenerate. 



The approach based on a numerical diagonal- 
ization of the space-discretized evolution operator 
IgI 17 1 effectively resolves all of the above issues. 
It allows the precise numerical calculation of a large 
number of energy eigenvalues and eigenstates of the 
single-particle Hamiltonian. Once calculated in the 
low-/3 regime, these data can be used to obtain the 
one-particle partition function for any value of the 
inverse temperature /3 in a numerically inexpensive 
way. Calculations based on this approach do not 
have restrictions in the high-/3 regime, where, in 
fact, they turn out to yield results with even bet- 
ter precision. In addition, the one-particle energy 
eigenfunctions obtained by the exact diagonaliza- 
tion will allow us to calculate local properties of 
Bose-Einstein condensates with very high accuracy. 

3. Numerical Calculation of Energy Eigen- 
values and Eigenstates 

The most efficient method for calculating prop- 
erties of few-body quantum systems is the di- 
rect diagonalization of the space-discretized prop- 
agator in imaginary time U{t) = exp{—tH/h) 
iSESbl, 22, 231. Here t represents the appro- 
priately chosen propagation time, for which we can 
calculate the matrix elements of the evolution oper- 
ator analytically within our path-integral approach 
in the = 1 approximation. 

On a given real-space grid defined by a set of 
coordinates xj = jA, where A denotes the spac- 
ing, j e [~L/ A, L/ A]'^ is a d-dimensional vector 
of integers, and L denotes the spatial cutoff, ma- 
trix elements of the propagator are just short-time 
transition amplitudes 





En/tuj±, k = ksEC 


n 


r = 


r = 1 





1.0009731351803 


0.1162667164134 


1 


2.0029165834022 


0.2674689968905 


2 


2.0029165834022 


0.2674689968905 


3 


3.0058275442161 


0.4426927375269 


4 


3.0058275442161 


0.4426927375270 


5 


3.0067964582067 


0.4725275724941 


6 


4.0097032385903 


0.6368178804983 


7 


4.0097032385903 


0.6368178804984 


8 


4.0116368851078 


0.6848142470356 


9 


4.0116368851078 


0.6848142470357 






En/flUj_, k = 10^ fcBEC 


n 


r = 


r ^ 1 





1.468486725893 


1.162667164134 


1 


3.213056378201 


2.674689968905 


2 


3.213056378201 


2.674689968905 


3 


5.163819069871 


4.426927375269 


4 


5.163819069871 


4.426927375270 


5 


5.406908088225 


4.725275724941 


6 


7.282930987460 


6.368178804982 


7 


7.282930987460 


6.368178804982 


8 


7.690584058915 


6.848142470357 


9 


7.690584058915 


6.848142470357 



[/ij(t) = (iA|J7(t)|jA), 



(15) 



Table 1: Lowest energy levels of the xy-pavt of the BEC 
potential JTJ for non-rotating (r = 0) and critically ro- 
tating (r = 1) condensate with the quartic anharmonicity 
k = ^BEC (top) and k = 10'' k^EC (bottom). They are 
obtained by using level p = 21 effective action with the dis- 
cretization parameters of Table[3] The spacing A was always 
chosen so that L/A = 100, and the propagation time was 
t = 0.2 for k = /cbec and t = 0.05 for k = 10^ /cbec- Errors 
are given by the precision of the last digit, typically 10~^^ to 
10-13, and are estimated by comparing the numerical results 
obtained with different discretization parameters. 



Note that throughout the paper we use dimension- 
less units, in which any energy is expressed in terms 
of huj±, while the length unit is the corresponding 
harmonic oscillator length y/h/Mjii)Uj±. 

The eigenvectors ?/'„(jA) of such a discretized 
evolution operator matrix correspond to the space- 
discretized eigenfunctions of the original Hamilto- 
nian, and the corresponding eigenvalues are re- 
lated to the eigenvalues of the single-particle Hamil- 
tonian by e"*^"/''. A detailed analysis of discretiza- 
tion errors of this method is given in Ref. 16]. It 
provides a practical algorithm for choosing optimal 
discretization parameters which minimizes the as- 
sociated errors. 



Another source of errors is related to the fact that 
for non-trivial potentials it is not possible to cal- 
culate exactly evolution operator matrix elements. 
Higher-order effective actions minimize the errors 
for a given propagation time t and allow an opti- 
mal choice of this parameter as well by taking into 
account other discretization errors, as was shown 
in Ref. [l3|. This approach is able to give very ac- 
curate energy eigenvalues even for moderate values 
of the propagation time t of the order 0.1. Tabled] 
presents the first several energy eigenvalues for the 
two-dimensional (xy-) part of the BEC potential ((1} 
for the non- rotating case (r = 0), as well as for the 





En/naj±, r = 1.04 


n 


fcBEC 


10' A:bec 





-0.6617041825660 


1.135693826206 


1 


-0.6465857464220 


2.628129903790 


2 


-0.6465857464220 


2.628129903790 


3 


-0.6032113415949 


4.363876633929 


4 


-0.6032113415948 


4.363876633929 


5 


-0.5349860004310 


4.667653582963 


6 


-0.5349860004309 


6.290444734007 


7 


-0.4451224795419 


6.290444734007 


8 


-0.4451224795419 


6.777210773169 


9 


-0.3362724309903 


6.777210773169 



Table 2: Lowest energy levels of the xj;-part of the BEC po- 
tential ((TJ for over-critically rotating (r = 1.04) condensate 
with the quartic anharmonicity k = /cbec s-^id k = 10^ ^bbc 
according to the same numerical procedure as in Tabled 



critically-rotating condensate (r = 1). The top ta- 
ble gives the energy spectrum of the potential with 
the anharmonicity k — fceEC used in the experiment 
while the bottom table shows the spectrum for 
the much larger anharmonicity k = 10'^ fcBEC- The 
degeneracies of numerically obtained eigenstates in 
all cases correspond to the expected structure of 
the spectrum, which can be deduced from the sym- 
metry of the problem. In addition to this, the in- 
teresting case of critical rotation [r — 1) allows 
a further verification of the numerical results. To 
this end we recall that the energy eigenvalues of a 
pure quartic oscillator, to which Vbec reduces in 
this case, are proportional to k^/'^ due to a spatial 
rescaling in the underlying Schrodinger equation. 
Therefore, we expect that the energy eigenvalues 
for k = 10' fcBEC are precisely 10 times larger than 
the corresponding eigenvalues for k = fceEC- Com- 
paring the rightmost columns in Table [1] we see ex- 
actly this scaling. This demonstrates conclusively 
that the presented method can be successfully ap- 
plied also in this deeply non-perturbative parame- 
ter regime. Furthermore, Table [5] gives the energy 
spectrum of an over-critically rotating (r = 1.04) 
condensate, illustrating that the same approach can 
be used in this delicate regime as well. 

With these results single-particle partition func- 
tions Zi{(3) can now be calculated according to 
Eq. ([7]). This is especially suitable for the low- 
temperature regime, when higher energy levels give 
a negligible contribution. Although the above de- 
scribed approach is able to accurately give several 
thousands of energy eigenvalues, their number is 
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60 80 

Figure 1: Numerically calculated density of states for xj/-part 
of the BEC potential with k = Aibec fo'' ^ critically rotating 
condensate, obtained by using level p = 21 effective action. 
The discretization parameters are L = 22.3, L/A = 100, and 
t = 0.2. The dashed line is the corresponding semiclassical 
approximation for the density of states. 



always necessarily limited. This is easily seen from 
Fig. [TJ where we compare the density of states for a 
critically rotating condensate with the correspond- 
ing semiclassical approximation for the density of 
states fl7'|. This comparison allows to estimate 
the maximal reliable two-dimensional energy eigen- 
value -Eniax which can be obtained numerically for 
a given set of discretization parameters. For ex- 
ample, from Fig. [1] we can estimate i^max ~ 90 for 
r = 1 with the discretization parameters L — 22.3, 
L/A = 100, and t = 0.2. Table [3] gives estimates 
for the maximal reliable energy eigenvalue for the 
anharmonicites k — fceEC and A: = 10' /cbec for 
several values of rotation frequencies. These re- 
sults are obtained from numerical calculations using 
the SPEEDUP codes [2J]. This table gives also an 
overview over those discretization parameters which 
were used for a numerical diagonalization of the 
BEC potential ([TJ in order to calculate both global 
and local properties of the condensate throughout 
the whole paper. 

In the low-temperature limit the finiteness of 
the number of known energy eigenstates does not 
present a problem. In fact, a precise knowledge 
of a large number of energy eigenvalues makes this 
approach a preferred method for a numerically ex- 
act treatment of low-temperature phenomena. On 
the other hand, the high-temperature regime, where 
thermal contributions of higher energy states play 
a significant role, is not treatable in the same way. 
This regime is usually not relevant for studies of 





k = /cbec 


k = 10^ A:bec 


r 




L 
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0.0 


140 


14.2 


190 


3.90 


0.2 


140 


14.4 


190 


3.90 


0.4 


140 


15.0 


180 


3.91 


0.6 


140 


16.3 


180 


3.92 


0.8 


130 


18.6 


180 


3.94 


1.0 


90 


22.3 


170 


3.96 


1.04 


90 


23.2 


170 


3.96 



Table 3: Maximal reliable numerically calculated energy 
eigenvalue Smax of the xi/-part of the BEC potential 1^ for 
different values of r = fl/uj±, estimated from comparing the 
numerically obtained density of states p{E) with the semi- 
classical approximation. The numerical diagonalization was 
done using level p = 21 effective action. The spacing A was 
always chosen so that L/A = 100, and the propagation time 
was t = 0.2 for k = Aibec ^^'^ * = 0.05 for k = 10^ fcaEC- 
The total number of reliable energy eigenstates is in all cases 
of the order of lO"*. 



BEC experiments, but we consider it for the sake 
of completeness. Furthermore, we want to demon- 
strate that the effective action approach can even be 
successfully used for studying this parameter range. 
In the case that the temperature is high enough, 
so that effects of higher energy eigenstates cannot 
be neglected, the inverse temperature /? becomes 
a small parameter. Thus, it becomes possible to 
calculate numerically the single-particle partition 
function as a sum of diagonal amplitudes, i.e. 

zm-T.U:i{m', (16) 

where A represents the spacing and, as before, the 
values of j are defined by j G [—L/A, L/A]'', with 
the spatial cutoff L chosen in such a way to ensure 
the localization of the evolution matrix within the 
interval [— L, L]*^. Note that this can be verified, 
since the evolution operator matrix elements are 
calculated using the analytic approximation (jl4p 
with level p effective potential 

where hl3 plays now the role of the imaginary time, 
and the influence of the kinetic term is not present 
since Ujj is a diagonal transition amplitude = 0). 



4. Finite Number of Energy Eigenvalues and 
Semiclassical Corrections 

In the previous section we have described a nu- 
merical approach which is capable of providing a 
large number of accurate energy eigenvalues for a 
general quantum system. Its major application is 
for few-body systems, where the complexity of the 
algorithm is very low. For instance, we are able 
to calculate typically lO'' energy eigenvalues for the 
considered BEC potential ((!]). In this section we 
discuss in more detail how the finiteness of nu- 
merically available energy eigenstates affects the 
calculation of thermodynamic properties of Bose- 
Einstein condensates. 

As outlined in Sec. [2j the information on single- 
particle eigenvalues is sufficient for calculating the 
condensation temperature according to Eq. (|10p . 
Below the condensation temperature, the ground- 
state occupancy follows from solving the equation 

oo 

iV = iVo + ^ [e'"'^^° Zi(m/3) - l] . (18) 

m—l 

In practical calculations, however, one is inevitably 
forced to restrict the sum over m in the cumulant 
expansion ((T8| to some finite cutoff M, resulting 
in the following approximation for the number of 
thermal atoms 

M oo 

N-No^J2Y. e-"«^"-^°) . (19) 

m— 1 n— 1 

Thus the ground-state occupancy A'o depends for 
each particle number N and temperature T also 
on this cutoff M. In particular, solving for 
the (inverse) condensation temperature, obtained 
by demanding iVo = 0, will yield Pc{M) with an 
explicit dependence on M . The exact condensation 
temperature /3c is only obtained in the limit M — > 
oo. 

Fig. [5] illustrates the M-dependence resulting 
from Eq. (1191) for both a non-rotating and a crit- 
ically rotating condensate. As expected, the sum 
saturates for high values of M to some finite num- 
ber N ~ Nq. By tuning the temperature in such a 
way that the sum saturates at the desired value of 
the total number of atoms N in the system, which 
implies A^o ~ 0, one is, in principle, able to extract 
the condensation temperature T^. 

Although the results in Fig. ^ suggest that this 
approach can be applied straightforwardly, a closer 
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Figure 2: Number of thermally excited atoms — A^o as a 
function of the cutoff M in the cumulant expansion II19I I . The 
results are given for a non-rotating condensate at T = 105.18 
nK and for a critically rotating condensate at T = 63.30 nK. 
The results are obtained by level p = 21 effective action, 
and all available numerical eigenstates are used to calculate 
N — Nq. The discretization parameters were L = 14.2 for 
r = and L = 22.3 for r = 1. In both cases the spacing was 
chosen according to L/A = 100, and the propagation time 
was t = 0.2. 



look at the results for numerically calculated values 
oi N — Nq reveals several problems that have to be 
addressed. At first we have to investigate how the 
results depend on the number of energy eigenstates 
used in the numerical calculation. Fig. [3] gives this 
dependence for a critically rotating condensate at 
its critical temperature Tc = 63.30 nK. We can see 
that the dependence on the maximal available two- 
dimensional energy eigenvalue E'max is quite signif- 
icant. The inset of this figure reveals another prob- 
lem: the value to which number Nq — N saturates 
depends in addition on the cumulant cutoff M, as 
explained earlier. While the M-dependence can be 
dealt with by using a very large value of the cu- 
mulant cutoff in numerical calculations, the depen- 
dence on the maximal energy eigenvalue £^max must 
be eliminated by taking into account a proper semi- 
classical correction to the single-particle partition 
functions. 

Namely, the finite number of energy eigenstates 
implies that the single-particle partition functions 
are only estimated by 



/3B„ 



(20) 
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Figure 3: Number of thermally excited atoms TV — A^o calcu- 
lated as a function of the maximal available two-dimensional 
energy eigenvalue Emax at T = 63.30 nK. The results are 
given for two different values of the cumulant cutoff M for a 
critically rotating condensate, with the same parameters as 
in Fig. [21 The horizontal line corresponds to the number of 
atoms N = 3 ■ 10^ in the experiment Q. 



A semiclassical correction to this value, can be cal- 
culated according to Ref. 10 as 

^|Pe-^«(x,p)e(iy(x,p)-iJ_), (21) 

where H{x,p) represents the classical Hamiltonian 
of the system, while Q denotes the Heaviside step- 
function. 

For the trap potential ([T|) we have in z-direction 
a pure harmonic potential which can be treated 
exactly. Therefore, we focus only on the two- 
dimensional problem in the xy-plane. In this case, 
the semiclassical correction for the single-particle 
partition function (|2ip can be expressed in terms of 
the complementary error function: 




where rimax corresponds to the value E'max of the 
numerically available maximal energy eigenvalue. 



xErfc 



When this semiclassical correction is taken into 
account, the numerical results show almost no de- 
pendence on i?max, as can be seen from Fig.U] Here 
we have used an excessively large value of the cu- 
mulant cutoff M — 10^ in order to completely elim- 
inate any M dependence. From the inset in this 




Figure 4: Number of thermally excited atoms N — No cal- 
culated as a function of -Emax with and without semiclas- 
sical corrections, calculated with any large cumulant cutoff 
M = 10* to eliminate the M-dependence. The results cor- 
respond to a critically rotating condensate with the same 
parameters as in Fig. [2] The horizontal line corresponds to 
A'^ = 301834 which represents the exact value at Tc = 63.30 
nK. 



Figure 5: Number of thermally excited atoms TV — A^o as 
a function of the temperature T for different values of the 
rotation frequency and the quartic anharmonicity k = /cbec ■ 
The discretization parameters are given in Table [S] and the 
results are calculated by taking into account semiclassical 
corrections. The dashed line corresponds to the number of 
atoms N = 3 ■ 10^ in the experiment [7|. For comparison, 
the full lines depict the semiclassical results from Ref. [lC|| . 



graph we also see that i?max must be chosen in ac- 
cordance with the value estimated in the previous 
section for the maximal reliable energy eigenvalue 
obtained by numerical diagonalization. If we use a 
value E'max larger than this, we will be underesti- 
mating the higher part of the energy spectra, and 
obtain incorrect results. For a critically rotating 
condensate with the anharmonicity k = A:bec the 
estimated value of i^max from Table [3] is around 90 
Sw_L, which agrees with the results from the inset 
of Fig. m If we use this value for i?max and calcu- 
late properties of the condensate using numerically 
obtained eigenstates below E'max with semiclassical 
corrections according to Eq. (P^ . we will obtain the 
exact results with very high accuracy. 

5. Global Properties of Rotating Bose- 
Einstein Condensates 

In this section we will apply this approach to cal- 
culate different global properties of rotating Bose- 
Einstein condensates. 

5.1. Condensation Temperature 

If we take into account semiclassical corrections 
as explained in the previous section, we can cal- 
culate, for instance, the condensation temperature 
of the condensate for different rotation frequencies. 
This implies that we have to find the temperature 



for which the number of thermal atoms N — Nq 
saturates precisely at the total number of atoms N . 
In practice, this works the other way around: for a 
given condensation temperature Tc we numerically 
calculate the particle number in the system using 
Eq. (jiyp . which gives the number of atoms in the 
system required for a condensation temperature to 
be equal to Tc. This procedure is implemented in 
Fig. [S] for several values of the rotation frequency ft 
in units of r = Q/uj±. For example, for Tc — 63.14 
nK we see that the corresponding number of parti- 
cles is iV = 3-10^, which coincides with the value for 
a critically rotating condensate in the experiment of 
Dalibard and co-workers Q- 

In principle, such a procedure is only applicable 
for low-accuracy calculations of the critical temper- 
ature, since otherwise one has to use very large val- 
ues of the cutoff M which would practically slow- 
down numerical calculations. If one is interested 
in more precise results, a suitable Af -dependence 
must be properly taken into account. In order to 
be able to efficiently extract the correct value of /3c , 
we will derive an analytical estimate of the asymp- 
totic error A/3c = Pc — l3c{M) which is introduced 
by the cutoff M. Note that always A/3c > 0, since 
/3c (Af) < /3c compensates the missing terms in the 
sum (fT9l) . 

If we insert /3c = /3c (M) + A/3c into Eq. (HH), the 
error A/3c can considered to be small for sufficiently 

9 



large value of the cutoff AI. By comparing Eq. 
with the exact expression (jlOp we obtain 



oo oo 

m=M+l n=l 
M 



A/3c J2 E"^(^»-^o)e" 



-m/3o(E„-£;o) 



(23) 



The term m(i?„ — i?o) within the sum can be ob- 
tained by setting Nq = and applying the partial 
derivative d/d^c to Eq. ([T^: 



1- A/3, 



m=M+l n=l 



-Bo) 



(24) 



Note that the derivative of the particle number N 
with respect to /3c is not equal to zero, since N is 
here effectively defined by the sum (jTU]). Therefore, 
we have instead 



dN 

Wc 



OO OO 

EE 

-m— 1 n— 1 



TO (£;„ - Eo) e-"'3'^(-^"--^") . (25) 



Clearly, the right-hand side is a negative quantity 
which does not depend on M . However, it does de- 
pend on /3c and the energy spectrum of the system. 

If the system is close to a d-dimensional harmonic 
oscillator, which is the case for the potential ([1]) 
with the small anharmonicity relevant for the ex- 
periment, for large values of M we have approxi- 
mately 



E E e~'"'^'^^^"~-^°^ 

m=M+l n=l 



1 _ g—fiahu: 



(26) 



where w denotes an effective harmonic frequency. 
For the case of a large anharmonicity, the effec- 
tive frequency w would depend on k, representing 
the harmonic expansion of the potential around its 
minimum. With such an estimate, Eq. ([M| reduces 
to 



A/3c 



dhuj 



(27) 
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The term {M + i) e-(*^+i)/3<:'*" in the denomina- 
tor of the second factor can be neglected for large 




Figure 6: Dependence of /3c on the cumulant cutoff M for 
an over-critically (r = 1.04) rotating condensate of N = 
3 ■ 10^ atoms of **^Rb with the quartic anharmonicity of the 
trap k = A:bec- The discretization parameters are given 
in Table |3] The dashed hne corresponds to a value of /3c 
obtained by fitting the numerical results to the function I I29I I. 
while the full line gives the fitted function f{M). 



enough values of the cutoff M, yielding as a simpli- 
fied version of the above expression: 



A/3c«- 



dN/dPc (1 - e-f^"'''^) 



(28) 



In order to use the derived estimates for A/3c, 
apparently one would already have to know the 
sought-after value of /3c as well as the difficult 
derivative dN/d/3c. However, in practical applica- 
tions this obstacle can be circumvented as follows. 
The expressions (P7|) and (^5]) can be used for fitting 
the numerical data for /3c(M) — /3c — A/3c, as is illus- 
trated in Fig. [6l In this standard approach, all un- 
known values are fit parameters, obtained numer- 
ically by the least-square method. Note that not 
only /3c is obtained by such a fitting procedure, but 
also other parameters, such as dN/d/3c, or the ef- 
fective harmonic frequency uj. The important point 
here is to capture the correct M-dependence, while 
all other parameters do not depend on it, so that 
they can be extracted by fitting. For example, in 
Fig. [5] we have used the fitting function 



/(M) = /3c 



cie 



-C2(M+1) 



1-HC3(M+ l)e-<=i(A'^+i) 



(29) 



which reproduces the numerical data quite accu- 
rately and gives high-precision results for the con- 
densation temperature Tc. The virtue of the de- 
rived estimates lies in the fact that they can be 
used to extract the information on the condensation 
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Figure 7: The condensation temperature as a function of 
the rotation frequency for the condensate of = 3 ■ 10^ and 
A'^ = 1-10* atoms of *^Rb, with the quartic anharmonicity of 
the trap k = /cbec- The discretization parameters are given 
in Table |3] The full lines correspond to the semiclassical 
approximation for Tc from Ref. [lOl . 



Figure 8: Ground-state occupancy Nq/N as a function of 
the temperature T for non-rotating and critically rotating 
condensate for different values of the total number of atoms. 
The quartic anharmonicity is fc = 10^ ^beCi ^^'^ the dis- 
cretization parameters are given in Table [3| The full lines 
depict the semiclassical results from Ref. |10|| . 



temperature even for moderate values of M, when 
a saturation is not yet achieved. This substantially 
speeds up the numerical calculation of condensation 
temperatures, especially when it has to be done for 
different values of potential parameters, such as the 
frequency ratio r = Q/uj^. 

Fig. [7] summarizes the numerical results for the 
condensation temperature Tc for the anharmonicity 
k = ksEC as well as the particle numbers = 3 - 10^ 
and A'' = 1 • 10^. If we compare the obtained numer- 
ical results with the semiclassical approximation 
from Ref. , we see that the agreement turns out 
to be relatively good for the undercritical regime, 
but it becomes worse for an overcritical rotation 
of the condensate. After presenting results for the 
ground-state occupancy, which were obtained from 
this approach in the next section, we will compare 
our numerically exact results with the semiclassi- 
cal approximation in more detail, and identify the 
parameter ranges where a full numerical treatment 
becomes necessary. 

5.2. Ground-State Occupancy 

The ground-state occupancy is the next impor- 
tant global property of Bose-Einstein condensates 
we will look into. Below the condensation tempera- 
ture a non-trivial fraction of atoms is in the ground 
state, thus yielding a macroscopic value of the oc- 
cupancy ratio Nq/N. 

Using the same approach as above, we can cal- 
culate the ground-state occupancy from Eq. ([T8|. 



After determining the ground-state energy Eq from 
an exact diagonalization of the evolution operator, 
we obtain 

AT 1 °° 

^ = l-iE[^'"'''°^i(-/^)-l] • (30) 

m— 1 

In order to calculate Nq/N, we need the full single- 
particle energy spectrum. For low temperatures, 
the large number of energy eigenstates obtained 
within the exact diagonalization is sufficient. In 
the sum of Eq. ([50]) we have again to introduce a 
cutoff M and to eliminate it by applying the meth- 
ods described in previous sections. To this end one 
uses either a very large value for the cutoff or one 
derives the appropriate finite correction term, and 
fits the results to the derived function. 

Fig. [5] presents numerical results for the ground- 
state occupancy of the condensate. The quartic 
anharmonicity of the trap is chosen to be fc = 
10'^ fceECj and the resuts are given for the nonrotat- 
ing case with the total number of atoms = 1 • 10^ 
and for critically rotating codensate with N = 5T0^ 
atoms. A comparison with the semiclassical results 
derived in Ref. shows that the deviations in- 
crease with larger temperatures and smaller parti- 
cle numbers. 

5.3. Comparison with Semiclassical Approximation 
Fig. ini depicts the errors of the semiclassically cal- 
culated condensation temperature, where the ex- 
act values are obtained by using the presented 
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Figure 9: Relative error of semiclassical results for the con- 
densation temperature [To| as a function of the rotation fre- 
quency fl in units of r = Q/u!±^ for TV = 3 ■ 10^ (top) 
and N = 1 ■ 10* (bottom). The quartic anharmonicity is 
k = fcBECi S'ld the discretization parameters are given in 
Table [3] The insets in both plots give the corresponding 
results for the large anharmonicity k = 10^ fcaEC- 



numerical approach. As we can see, the agree- 
ment is relatively good for large particle numbers 
and small anharmonicity if the condensate rotates 
under-critically. The error in this case is of the or- 
der of 1 % to 1.5 %, and turns out to be minimal for 
a critical rotation. However, the error significantly 
increases for an overcritical rotation up to almost 
3.5 % for r — 1.1. Therefore, while the semiclassical 
approximation is acceptable for undercritical rota- 
tion, in the overcritical regime a numerical treat- 
ment becomes necessary. This is even more pro- 
nounced if we decrease the particle number to 10^, 
which is quite typical for many BEC experiments. 
In that case, semiclassical results already have an 
error of the order of 20 %. For large anharmonicity 
the rotation effect is not so important, as we can 
see from insets on both graphs in Fig. |9l However, 



increasing the particle number to 10^ makes a nu- 
merical treatment indispensable, since the errors of 
the semiclassical results amount up to 5 %. 



6. Local Properties of 
Einstein Condensates 



Rotating Bose- 



Local properties of ultra-cold quantum gases 
are ubiquitously used to observe and study the 
phenomenon of Bose-Einstein condensation. The 
prominent peak in time-of-fiight absorption pic- 
tures, which appears suddenly when the temper- 
ature is decreased below Tc, is a clear signature for 
the occurrence of a BEC phase transition. It is ex- 
perimentally used to measure the thermodynamic 
properties of the condensate. In this section we will 
show how the presented numerical approach can be 
applied to calculate both the density profiles and 
the time-of-fiight absorption imaging profiles. 

6.1. Density Profiles 

The two-point propagator p(xi,X2) = 
(^^(xi)*I'(x2)) defines via its diagonal element, i.e. 
n(x) = /9(x, x), the density profile of atoms in a 
trap. For the ideal Bose gas, the density profile 
can be written as 



i(x) = 7Vo|V'o(x)P + ^iV„|Vn(x)| 



(31) 



where the second term represents the thermal con- 
tribution to the density profile. Furthermore, the 
quantities ^„ represent single-particle eigenstates, 
while the occupancies iV„ with n > 1 are given by 
the Bose-Einstein distribution 



1 



c/3(E„-Bo 



1 



(32) 



Having at our disposal numerically calculated en- 
ergy eigenvalues and eigenfunctions, we can calcu- 
late the density profile of the condensate. In or- 
der to do so, we first have to obtain the ground- 
state occupancy number iVo using the approach de- 
scribed in the previous section. Once this is done, 
Eq. (j3ip allows to calculate the density profile. 
In view of a comparison with absorption imaging, 
which always produces two-dimensional profiles, we 
have to integrate our numerically determined three- 
dimensional particle density ?t.(x) along the imag- 
ing axis. Fig. [10] presents typical results for the 
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Figure 10: Density profile in xy-plane for a non-rotating 
(top) and a critically rotating (bottom) condensate of N = 
3 ■ 10^ atoms of *^Rb with the anharmonicity k = fcBEC at 
T = 30 nK. The dimensionless unit length on both graphs 
corresponds to 1.34 fim, i.e. the linear size of the profile is 
approximately 16.1 fim (top) and 32.2 fim (bottom). The 
discretization parameters are given in Table [3] 



resulting density profiles of Bosc-Einstcin conden- 
sates for both the non-rotating and the critically- 
rotating case. Obviously, a rotation of the conden- 
sate leads to an effective spreading due to the ap- 
pearance of a centrifugal potential. 

Although this approach is sufficient for treating 
the low-temperature regime, where the condensate 
is present, we emphasize that the same method can 
also be used to deal with the thermal regime, when 
the temperature is increased beyond T^. For even 
higher temperatures, when the number of energy 
eigenstates, that need to be taken into account, ex- 
ceeds the number of numerically accessible eigen- 
states, the presented approach can be extended in 
a similar way as the partition function was calcu- 
lated previously as a sum of diagonal transition 
amplitudes. Using the cumulant expansion of oc- 
cupancies and the spectral decomposition of ther- 
mal transition amplitudes, the density profile can 
be written for high enough temperatures as 



n(x) =7Vo|V^o(x)|2 

+ [e'"''^M(x, 0; X, m/3n) 



Here A(x, 0; x, mj3h) represents the imaginary-time 
amplitude for a single-particle transition from the 
position X in the initial imaginary time t = Q io the 
position X in the final imaginary time t = mj3h. 

While both definitions (I5T]) and p3l) are math- 
ematically equivalent in the case when one is able 
to calculate infinitely many energy eigenstates and 
amplitudes for an arbitrary propagation time, the 
first one is more suitable for low temperatures, 
when the number of relevant energy eigenstates 
is moderate, and the second one is suitable for 
high temperatures, when the imaginary propaga- 
tion time hj3 is small, and the short-time expansion 
can be successfully applied. 

6.2. Time- of- Flight Graphs for BECs 

In typical BEG experiments, a trapping potential 
is switched off and the gas is allowed to expand 
freely during a short flight time t which is of the 
order of several tens of miliseconds. Afterwards an 
absorption picture is taken which maps the density 
profile to the plane perpendicular to the laser beam. 
For the ideal Bose condensate, the density profile 
after time t is given by 

n(x,t) Aro|V'o(x,t)P + ^A^„|^„(x,t)|2, (34) 

r!>l 

where the density profile has to be integrated along 
the imaging axis, and the eigenstates ?A„(x, are 
propagated according to the free Hamiltonian, con- 
taining only the kinetic term, since the trapping po- 
tential is switched off. If the energy eigenstates are 
available exactly, either analytically or numerically, 
their propagation time can be calculated by per- 
forming two consecutive Fourier transformations: 



-0n(x,t) = 



(27r)3 



gZ[k.(r-R)-^ktl^^^(R)^ (35) 



|Vo(x)n .(33) 



where the term e"'"''* accounts for a free-particle 
propagation in k-space. In practical applications, 
when the energy eigenstates are calculated by a nu- 
merical diagonalization of space-discretized transi- 
tion amplitudes, the natural way to calculate the 
above free-particle time evolution is to use Fast 
Fourier Transform (FFT) numerical libraries. 

For high temperatures we can use a mathemat- 
ically equivalent definition of the density profile 
which is derived again from using the cumulant ex- 
pansion of occupancy numbers and the spectral de- 
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Figure 11: Time-of-flight absorption density profiles in xy- 
plane for an over-critically rotating (r = 1.04) condensate 
of Af = 3 • 10^ atoms of **^Rb with the anharmonicity k = 
fcBEC at T = 30 nK. The flight time, designated as TOF, 
is given at each plot. The dimensionless unit length on all 
graphs corresponds to 1.34 fim and the linear size of profiles 
is approximately 53.6 ^m. The discretization parameters are 
given in Table |3] 



composition of transition amplitudes: 

^m^Eo f d^^l ^'k2 d^Xi d^X2 



E 

m>l 



(2n) 



„i[(ki-k2)-x-ki-Xi+k2-X2-(cJki -i^kj )*] 



xA(Xi,0;X2,m/3/i) - |V'o(x,i)l 



(36) 



In both approaches it is first necessary to calcu- 
late the ground-state energy Eq and the eigenfunc- 
tion 'i/'o(x), as well as the ground-state occupancy 
A^o- If we rely on Eqs. and ([55)1 to calculate 
time-of-flight graphs, we have to calculate as many 
eigenstates as possible by numerical diagonaliza- 
tion. Conversely, if it is possible to use directly 
Eq. ((36)) . we can apply the effective action short- 
time expansion of thermal transition amplitudes. 
In both cases FFT is ideally suited for calculating 
time-of-flight graphs. 

6.3. Overcritical Rotation 

The case of critical and overcritical rotation r > 1 
is realized in the Paris experiment by introducing 
the anharmonic part of the potential ([1]), so that the 
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Figure 12: Condensate density at the origin of xj;-plane as 
a function of the time of flight (TOF) for the condensate of 
N = 3-10^ atoms of *''Rb at T = 30 nK for several rotation 
frequencies fl in units of r = Q/lu±. The quartic anhar- 
monicity is fe = fcBEC and the discretization parameters are 
given in Table [3] 



condensate is confined even when the harmonic part 
of the trapping potential is completely compensated 
or overcompensated by the rotation. The experi- 
mental realization of this delicate balance was diffi- 
cult to achieve, but nevertheless when the conden- 
sate was successfully confined while rotating over- 
critically, the measurements of its properties can 
be done using the standard techniques, including 
absorption imaging. Within the semiclassical ap- 
proach one has to carefully consider this situation, 
since the chemical potential is defined by the mini- 
mum of the potential, and now cannot be simply set 
to zero anymore 10 . In our numerical approach, 
however, the implementation of the methods de- 
scribed in previous sections is straightforward even 
for overcritical rotation. First one calculates en- 
ergy eigenvalues and eigenstates using exact diago- 
nalization, yielding negative values for the first sev- 
eral eigenstates. Table[l]shows the resulting energy 
spectrum of an over-critically rotating condensate 
(r = 1.04) for the experimental value of the an- 
harmonicity /cbeCj as well as for the case of large 
anharmonicity 10^ ^bec- 

Condensation temperature and other global 
properties as well as local properties of overcriti- 
cally rotating condensates can also be calculated as 
before. Fig. [TT] gives the time-of-fiight absorption 
imaging sequence in the xy-plane for an overcriti- 
cally (r = 1.04) rotating Bose-Einstein condensate 
with the anharmonicity k — /cbec and the particle 
number = 3 • lO'^ at T = 30 nK, exhibiting an 
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interesting behavior. The initial density profile has 
a minimum at the origin, due to the shape of the 
anharmonic potential. The free expansion of the 
condensate leads to an increase in the particle den- 
sity at the origin, and only afterwards the conden- 
sate density profile expands monotonically. Fig. [T^] 
presents the time dependence of the particle den- 
sity at the origin for varying rotation frequencies, 
parametrized by the ratio r = ^l/uj^. We read off 
that approaching the critical rotation slows down 
the expansion of the condensate. For overcritical 
rotation this is even more pronounced, due to the 
appearance of the peak in the particle density for 
the expansion time t > Q. This leads to an expan- 
sion which is typically an order of magnitude slower 
for the rotation with r > 1. 

7. Conclusions 

In this Letter a new method for numerically cal- 
culating short-time transition amplitudes based on 
the effective action approach is applied to the study 
of ideal Bose gases. Earlier derived higher-order dis- 
cretized effective actions are used for an efficient nu- 
merical calculation of both global and local proper- 
ties of fast-rotating Bose-Einstein condensates. To 
this end we have calculated large numbers of single- 
particle eigenvalues and eigenstates using an exact 
numerical diagonalization of the space-discretized 
evolution operator matrix. Using this information, 
we have calculated the condensation temperature 
and the ground-state occupancy of the condensate, 
as well as density profiles and time-of-flight absorp- 
tion graphs. We have also shown that a critical and 
an overcritical rotation can be studied using the 
presented numerical approach, and that it leads to 
a substantial increase in the time scale for the free 
expansion of the condensate after the trapping po- 
tential is switched off. Finally, we note that our 
approach can also be used for numerical studies of 
properties of rotating ultra-cold Fermi gases pH ). 
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